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ABSTRACT 

We  have  compiled  an  extensive  data  set  on  potential  parts  of  the  Monoceros  tidal  stream  and  performed  an 
exhaustive  survey  of  dwarf  galaxy  semianalytic  orbits  in  order  to  constrain  its  orbital  properties.  The  best-fit  orbits  are 
subsequently  realized  as  self-consistent  A-body  simulations  in  order  to  reproduce  the  spatial  and  velocity  distribution 
of  satellite  debris.  We  find  that  all  kinematic  and  geometric  constraints  can  be  fit  by  a  single  stream  allowing  for 
multiple  wraps.  The  orbital  eccentricity  and  inclination  of  the  progenitor  are  strongly  constrained  to  be  e  =  0.10  zb 
0.05  and  i  =  25°  d=  5°.  Ten  new  estimates  of  proper  motions  from  the  Sloan  Digital  Sky  Survey  clearly  exclude  all 
retrograde  orbits.  Particles  lost  by  the  satellite  populate  two  nearly  concentric  rings,  naturally  explaining  the  detection 
of  stream  stars  at  both  6-8  kpc  (Ibata  et  al.;  Newberg  et  al.)  and  12-18  kpc  (the  Tri/And  stream;  Rocha-Pinto  et  al.) 
from  the  Sun.  We  have  attempted  to  predict  the  present  location  of  the  Monoceros  stream  progenitor  using  different 
information:  (1)  the  kinematical  and  spatial  distribution  of  detections,  and  (2)  the  different  mean  metallicity  in  the 
inner  and  the  outer  rings.  Because  of  the  lack  of  observational  data  in  the  whole  range  of  Galactic  latitudes,  the 
geometrical/kinematical  constraints  lead  to  a  wide  range  of  possible  locations.  By  associating  older  parts  of  the  model 
stream  with  lower  metallicity  parts  of  the  observed  data,  we  argue  in  favor  of  a  current  location  of  /  ~  245°, 
b  ~  —  1 8°,  with  a  distance  to  the  Sun  rs  ~  1 5  kpc.  The  mass  of  the  progenitor  has  been  poorly  constrained  because  of 
the  slow  orbital  decay.  Similar  fits  have  been  obtained  for  masses  (3-9)  x  108  M0.  We  have  analyzed  the  possible 
common  origin  of  the  Canis  Major  dwarf  and  the  Monoceros  stream.  The  Canis  Major  dwarf  moves  on  a  pro  grade, 
nearly  circular  orbit  (e  ~  0. 1 6)  in  the  Milky  Way  disk  (i  ~  4^4  deg).  This  orbital  inclination  is  too  low  to  account  for 
the  large  vertical  dispersion  of  stream  stars.  However,  the  bimodal  distribution  of  radial  velocities  in  the  central  region 
found  by  Martin  et  al.  probably  indicates  that  their  selection  criteria  for  identifying  dwarf  stars  lead  to  a  contamination 
of  background  stars.  In  that  case,  the  kinematical  data  outlined  above  might  result  in  an  underestimate  of  the  orbital 
inclination.  Finally,  the  distance  estimation  to  Canis  Major  dwarf  is  around  a  factor  of  2  smaller  than  that  obtained 
from  our  model.  Unfortunately,  the  possible  identification  of  the  Monoceros  stream  progenitor  in  Canis  Major 
remains  unclear. 

Subject  headings:  galaxies:  evolution  —  galaxies:  formation  —  galaxies:  halos  —  galaxies:  structure 
Online  material:  color  figures 


1.  INTRODUCTION 

In  a  ACDM  universe,  the  inner  regions  of  massive  galaxies 
like  the  Milky  Way  gain  a  large  fraction  of  their  mass  through 
tidal  disruption  and  accretion  of  a  large  number  of  low-mass 
fragments  (e.g.,  van  den  Bosch  et  al.  2004).  The  fossil  records  of 
these  merging  processes  may  be  observable  nowadays  in  the 
form  of  long  tidal  streams  or  large-scale  stellar  substructures 
around  the  parent  spiral  galaxies.  Numerical  galaxy  formation 
simulations  show  that  dynamical  friction  brings  substructures 
from  outer  halo  regions  to  the  neighborhood  of  the  parent  gal¬ 
axies’  disks.  If  halos  are  flattened  (oblate)  and  their  axisymmetry 
plane  is  that  of  the  disk,  orbits  of  nonpolar  satellite  galaxies  tend 
to  become  coplanar  and  circularize  with  time  (Penarrubia  et  al. 
2004).  In  this  scenario,  the  tidal  debris  of  several  disrupted  sat¬ 
ellite  galaxies  might  have  contributed  to  the  formation  of  the 
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stellar  disk  if  they  followed  nearly  circular  orbits  with  a  low 
orbital  inclination  at  late  times  of  their  evolution  (Navarro  2004 
and  references  therein). 

The  Milky  Way  is  an  important  laboratory  in  which  the  pre¬ 
dictions  of  this  cosmological  scenario  can  be  tested.  In  the  last 
decade,  large-scale  surveys  have  proved  the  existence  of  tidal 
streams  (Sagittarius  [Ibata  et  al.  1994]  and  Monoceros  [Newberg 
et  al.  2002])  in  our  Galaxy,  providing  strong  observational  evi¬ 
dence  that  disruption  of  dwarf  satellites  contributes  to  the  as¬ 
sembly  of  some  components  of  our  Galaxy.  These  tidal  streams 
offer  a  unique  opportunity  to  study  accretion  events  in  consider¬ 
able  detail  using  the  chemical,  kinematic,  and  spatial  distribution 
of  tidal  stream’s  stars,  which  can  be  directly  compared  against 
predictions  of  A-body  simulations  of  merging  events  (e.g.,  Law 
et  al.  2005). 

Recently,  the  Sloan  Digitized  Sky  Survey  (SDSS)  team  reported 
the  discovery  of  a  coherent  ringlike  structure  at  low  Galactic 
latitude  spanning  about  100°  in  the  sky  (Newberg  et  al.  2002; 
Yanny  et  al.  2003).  Follow-up  observations  (Ibata  et  al.  2003) 
found  that  this  structure  of  low-metallicity  stars  surrounds  the 
Galactic  disk  at  Galactocentric  distances  from  ~15  to  ~20  kpc. 
Different  scenarios  have  been  proposed  to  explain  the  nature 
of  this  feature  in  the  outer  Galactic  disk  (see  Helmi  et  al.  2003): 

(1)  a  new  tidal  stream  (Yanny  et  al.  2003;  Ibata  et  al.  2003),  and 

(2)  a  stellar  ring  resulting  from  perturbations  in  the  disk  similar 
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TABLE  1 

Observational  Constraints 


Authors 

l-b  Range 
(deg) 

Type 

vra 

nh 

Nc 

Symbol 

Yanny  et  al.  (2003) . 

[182,  225],  [+28,  -27] 

CMD 

Yes 

No 

4 

Filled  triangle 

Ibata  et  al.  (2003) . 

[122,  218],  [+30,-25] 

CMDd 

No 

No 

14 

Open  triangle 

Crane  et  al.  (2003) . 

[157,  242],  [  +  38,  -15] 

M  giant 

Yes 

Yese 

58 

Filled  triangle 

H.  J.  Newberg  et  al.  (2005,  in  preparation) . 

[110,  225],  [  +  37,  -32] 

CMDd 

No 

No 

22 

Open  square 

Rocha-Pinto  et  al.  (2003) . 

[117,  157],  [  +  38,  -25] 

M  giantf 

Yes 

No 

31 

Filled  square 

a  Radial  velocity  measurements. 
b  Proper-motion  measurements. 
c  Number  of  detections. 

d  Possible  confusion  with  thick-disk  stars  (see  §  5.4). 
e  Only  for  10  stars. 
f  Tri/And  stream. 


to  the  remnants  of  ancient  warps  (Ibata  et  al.  2003).  Tracing  this 
structure  with  Two  Micron  All  Sky  Survey  (2MASS)  M-giant  stars, 
Rocha-Pinto  et  al.  (2003)  concluded  that  its  structural  charac¬ 
teristics  are  consistent  with  an  interpretation  of  this  structure 
as  the  fossil  of  a  merging  dwarf  galaxy,  similar  to  the  Sagittarius 
dwarf  galaxy  (Ibata  et  al.  1994),  in  the  Galactic  plane.  Interest¬ 
ingly,  Frinchaboy  et  al.  (2004)  and  Crane  et  al.  (2003)  also  sug¬ 
gested  some  nearby  Galactic  open  and  globular  clusters  with 
coordinated  heliocentric  radial  velocities,  indicating  a  possible 
common  origin  with  the  tidal  stream. 

Unlike  the  Sagittarius  dwarf  galaxy,  the  Monoceros  stream 
has  been  detected  prior  to  locating  the  main  body  of  the  parent 
galaxy.  The  location  is  still  controversial.  From  the  analysis 
of  the  2MASS  catalog,  Martin  et  al.  (2004a)  reported  a  strong 
elliptical-shaped  stellar  overdensity  in  the  constellation  of  Canis 
Major,  which  is  coincident  in  position  and  radial  velocity  with  a 
small  group  of  four  globular  clusters  (see  §  5.3).  Using  a  fairly 
simple  model,  they  identify  the  CMa  overdensity  as  the  main 
body  of  the  progenitor  dwarf  galaxy  of  the  Monoceros  stream 
(named  the  CMa  dwarf).  Bellazzini  et  al.  (2004)  presented  color- 
magnitude  diagrams  in  the  surroundings  of  the  CMa  location, 
concluding  that  the  system  is  situated  at  8  it  1  kpc  from  the  Sun 
and  that  it  is  composed  of  a  metal-rich,  intermediate-age  popu¬ 
lation.  Momany  et  al.  (2004),  however,  comparing  proper  mo¬ 
tions  and  radial  velocities  of  CMa  M-giant  stars,  obtain  that  the 
overdensity  in  this  region  mimics  the  thick-disk  kinematics. 
Moreover,  they  find  that  the  star  counts  in  that  region  are  those 
expected  in  the  standard  Milky  Way  model  if  one  takes  into 
account  the  warp  and  flare  of  the  disk.  However,  a  deep  color- 
magnitude  diagram  of  the  center  of  the  CMa  overdensity  by 
Martinez-Delgado  et  al.  (2004)  shows  a  very  well-defined  main 
sequence  consistent  with  a  limited  extent  in  distance,  confirming 
that  this  stellar  population  is  associated  with  a  distinct,  possibly 
still  bound  stellar  system  with  properties  (surface  brightness,  ab¬ 
solute  magnitude,  stellar  content)  compatible  with  those  of  Local 
group  dwarf  spheroidal  galaxies.  On  the  basis  of  these  last  re¬ 
sults,  the  CMa  overdensity  will  be  referred  to  the  CMa  dwarf 
galaxy  in  this  work. 

In  this  paper  we  present  the  results  of  an  extensive  search 
of  possible  disruption  scenarios  of  satellite  galaxies  to  explore 
(1)  how  many  of  the  “overdensity  signatures”  can  be  attributed 
to  a  single  stream  and  (2)  to  constrain  its  progenitor  orbit.  With 
an  objective  criterion  we  select  those  orbits  that  reproduce  the  spa¬ 
tial  and  kinematical  distribution  of  candidate  debris.  Beyond  explor¬ 
ing  the  possible  progenitor  location,  the  results  will  (1)  help  to 
carry  out  new  surveys  in  different  regions  of  the  sky  and  (2)  con¬ 


strain  the  thick-disk  formation  history  through  the  distribution 
evolution  of  stripped  stars. 

This  paper  has  been  organized  as  follows.  In  §  2  we  compile 
the  observational  data  available  so  far.  Section  3  outlines  the  gal¬ 
axy  and  satellite  models  that  we  use  in  our  investigation.  We  also 
describe  our  semianalytic  orbit  algorithm  and  the  results  of  the 
orbital  parameter  survey.  We  obtain  the  distribution  of  debris  via 
the  TV-body  algorithm  sketched  in  §  3.4.  In  §  4  we  compare  the 
best-fit  orbit  with  the  observational  data,  and  in  §  5  we  comment 
on  different  implications  of  our  results. 

2.  OBSERVATIONAL  DATA 

For  the  present  work,  we  have  attempted  a  comprehensive 
compilation  of  observational  constraints  on  the  geometry  and  ki¬ 
nematics  of  the  Monoceros  stream.  Given  the  great  deal  of  activity 
directed  toward  detecting  tidal  streams  in  general,  the  diversity 
of  available  data  sets,  and  in  particular  the  great  interest  in  the 
Monoceros  stream  spawned  by  the  original  discovery,  the  con¬ 
straints  are  both  numerous  and  inhomogeneous. 

Typically,  detections  of  overdensities  have  resulted  from  pho¬ 
tometry,  finding  an  excess  of  main-sequence  stars  or  M  giants  at 
a  given  apparent  magnitude.  Depending  on  the  tracer  population, 
this  approach  provides  distance  estimates  of  differing  quality. 
Radial  velocities  are  available  for  a  subset  of  stars  within  a  subset 
of  directions  with  detected  overdensities.  In  principle,  the  chem¬ 
ical  composition  of  stars  can  also  serve  as  a  constraint,  as  it  can 
be  expected  to  vary  continuously  along  the  stream,  presumably 
with  decreasing  metallicity  in  the  most  loosely  bound  (and  hence 
first  lost)  material.  Finally,  proper  motions  (which  typically  can 
be  measured  to  A /jl  ~  3-4  mas  yr-1)  can  provide  some  rough 
constraints;  their  precision  in  physical  units  is,  however,  only 
At;perp  =  A fixrs  x4.74,  where  rs  is  the  distance  from  the  Sun 
and  4.74  is  the  factor  that  converts  (kpc  mas  yr-1)  into  (km  s-1). 

We  have  compiled  both  published  constraints  and  those  avail¬ 
able  to  us  but  still  in  the  process  of  publication.  A  few  constraints 
were  derived  specifically  for  the  present  paper.  We  have  sum¬ 
marized  the  constraints,  their  nature  (photometric,  kinematic, 
etc.),  and  their  sources  in  Table  1.  In  the  following  plots  we  use 
different  symbols  to  distinguish  between  different  data  sources. 
In  order  to  simplify  our  figures  we  have  made  an  exception  for 
those  data  with  available  radial  velocities,  which  we  have  been 
plotted  everywhere  with  filled  triangles. 

In  addition,  we  have  proper  motions  of  confirmed  stream 
star  members  selected  from  the  radial  velocity  compilation  by 
Crane  et  al.  (2003).  Proper  motions  are  determined  by  combining 
recalibrated  USNO-Bl.O  positions  (Monet  et  al.  2003)  with 
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TABLE  2 

Proper  Motions  of  Monoceros  Stream  Stars 


/ 

(deg) 

b 

(deg) 

rs 

(kpc) 

vr 

(km  s-1) 

ni 

(mas  yr-1) 

n-b 

(mas  yr-1) 

186.382 . 

23.910 

11.2 

-6.8  ±  3.9 

2.6 

-1.9 

189.316 . 

23.251 

12.6 

149.4  ±  5.6 

0.6 

-4.4 

186.894 . 

24.181 

10.0 

5.4  ±  3.7 

-5.5 

-13.9 

189.741 . 

23.345 

12.7 

0.2  ±  6.0 

6.0 

-1.0 

198.778 . 

25.063 

11.9 

49.1  ±  2.5 

-1.2 

2.6 

224.108 . 

21.344 

10.9 

82.8  ±  3.4 

1.0 

-6.3 

178.371 . 

36.786 

10.4 

-19.7  ±  2.0 

6.1 

1.5 

186.992 . 

38.816 

12.3 

42.3  ±  2.1 

0.4 

3.0 

221.989 . 

29.900 

12.3 

86.1  ±  3.2 

-3.3 

-6.2 

223.126 . 

32.262 

11.8 

55.1  ±  2.7 

-3.5 

1.9 

Notes. — Positions  in  Galactic  coordinates  of  the  Monoceros  stream  stars  with 
measured  proper  motions.  The  errors  in  the  heliocentric  distance  rs  are  about  25% 
of  the  value,  whereas  for  /i/  and  gb  they  have  been  estimated  to  be  3.5  mas  yr-1. 


nally,  we  present  the  method  used  to  constrain  the  orbital  prop¬ 
erties  of  a  possible  Monoceros  tidal  stream  progenitor,  which 
can  be  divided  in  two  steps.  (1)  Satellite  orbits  are  calculated  us¬ 
ing  a  well-tested  semianalytic  algorithm  in  order  to  perform  a 
survey  of  our  large  parameter  space.  Subsequently,  the  orbit  col¬ 
lection  has  been  compared  against  the  observational  sample  to 
determine  the  parameter  sets  that  provide  the  best  fits.  (2)  Once 
the  orbital  constraints  are  known,  we  carry  out  TV-body  simula¬ 
tions  to  analyze  the  spatial  and  kinematical  distribution  of  debris. 

3.1.  Galaxy  and  Dwarf  Satellite  Models 

Our  dynamical  model  for  the  Milky  Way  follows  Hemquist 
(1993),  where  (1)  the  disk  is  exponential  in  the  radial  direction 
and  isothermal  in  the  vertical  direction;  (2)  the  bulge  is  described 
by  a  Hemquist  model  (Hemquist  1990)  with  spherical  symme¬ 
try;  and  (3)  the  halo  follows  a  flattened,  nonsingular  isothermal 
profile  with  given  core  and  cutoff  radius. 

For  the  density  distributions  of  the  disk  we  take 


those  from  the  SDSS7  catalog,  as  detailed  in  Munn  et  al.  (2004). 
The  proper  motions  of  our  target  stars  are  given  in  Table  2.  It  is 
necessary  to  clarify  that  the  measurements  of  proper  motions 
have  not  been  included  in  the  orbital  fit  and  will  be  used  only  to 
discriminate  the  sense  of  motion  of  the  progenitor’s  orbit. 

In  addition,  we  are  faced  with  the  problem  of  which  con¬ 
straints  (and  on  which  candidate  sections  of  the  Monoceros 
stream)  to  include  in  our  modeling.  Guided  by  the  goal  of  finding 
the  largest  number  of  stream  portions  that  may  have  arisen  from 
a  single  dismption  event,  we  have  used  an  iterative  procedure. 
Starting  with  the  original  stream  detections  (Newberg  et  al.  2002; 
Yanny  et  al.  2003)  and  the  comprehensive  kinematic  survey  of 
M  stars  (Crane  et  al.  2003)  spanning  100°  in  the  sky.  Initial  mod¬ 
eling  of  these  constraints  made  it  clear  that  other,  independently 
found  overdensities  are  most  likely  also  part  of  the  stream.  In 
particular,  the  Tri/And  stream  (Majewski  et  al.  2004),  a  more  dis¬ 
tant  metal-poor  stellar  stream,  showed  a  location  in  the  sky  and 
a  radial  velocity  curve  fairly  similar  to  the  predictions  of  our 
first-iteration  model.  Therefore,  the  available  data  (including 
radial  velocities  from  Rocha-Pinto  et  al.  2003)  on  the  Tri/And 
stream  were  then  included  as  inputs  in  a  second  iteration  of  our 
survey  of  the  best-candidate  orbit  to  better  constrain  its  prop¬ 
erties  and  reduce  the  number  of  possible  scenarios. 

Finally,  we  explicitly  test  whether  the  CMa  dwarf  is  likely  part 
of  the  stream,  presumably  the  “parent”  of  the  tidal  debris.  We  do 
this  by  omitting  CMa  in  the  first  modeling  and  then  comparing 
its  orbit  with  that  of  our  model  for  a  possible  Monoceros  stream 
progenitor. 

3.  FITTING  NUMERIC AF  MODEFS  TO 
OB  SERVATIONAF  DATA 

In  this  section  we  outline  how  we  perform  the  orbital  analysis 
of  dwarf  satellites  and  the  method  used  to  fit  numerical  orbits  to 
the  observational  data. 

We  first  describe  a  Milky  Way  model  that  matches  the  mass 
distribution  of  our  Galaxy.  Subsequently,  we  discuss  the  param¬ 
eter  space  that  we  must  cover  in  order  to  determine  the  orbital 
characteristics  of  a  possible  Monoceros  stream  progenitor.  Fi- 


Pd(R,z )  =  MJ2  exp  (-R/Rj)sech2(z/z0),  (1) 

^Rdz o 

where  Mg  =  5.60  x  1010  M0  is  the  disk  mass,  z0  =  0.70  kpc 
is  the  vertical  thickness,  and  Rg  =  3.50  kpc  is  the  exponential 
scale  length  in  the  radial  direction.  The  mass  profile  decays  ex¬ 
ponentially  with  R  and  is  composed  of  isothermal  sheets  along 
the  vertical  direction.  Velocities  are  assumed  to  have  a  Gaussian 
distribution. 

For  the  bulge  we  adopt  the  spherical  Hemquist  profile 
(Hemquist  1990), 


a 

27 r  r{r  +  af  ’ 


(2) 


where  =  1 .86  x  1010  M0  is  the  bulge  mass  and  a  =  0.53  kpc 
is  the  spherical  scale  length.  This  analytical  profile  fits  the  de 
Vaucouleurs  law  (de  Vaucouleurs  1948).  The  velocity  field  is 
constmcted  from  the  Jeans  equations  by  assuming  isotropic 
Gaussian  velocity  distributions  at  each  radial  distance  (Hemquist 
1993). 

We  use  a  nonsingular  isothermal  profile  for  the  dark  matter 
halo  (DMH), 


Mh<x  exp(-m2/r2ut) 
2i r3/2rcut  m2 -\-  7  2 


(3) 


where 


m2=R2+z2/ql  (4) 

in  cylindrical  coordinates;  is  the  halo  density  flattening,  M \  = 
7.84  x  1011  Mq  is  the  DMH  mass,  rcut  =  84.00  kpc  the  cutoff 
radius,  7  =  3.50  kpc  the  core  radius,  and 

a  ={\  -  V^r/^ex p(/J2)[l  -  erf(/J)]}-1 

=  1  +  ~\~  (7r  —  2)/J2  +  0(/73),  (5) 


7  SDSS  (York  et  al.  2000)  is  an  imaging  and  spectroscopic  survey  that  will 
eventually  cover  of  the  sky.  Drift-scan  imaging  in  the  five  SDSS  bandpasses 
(u,  g,  r,  i,  z )  (Fukugita  et  al.  1996;  Gunn  et  al.  1998;  Hogg  et  al.  2001)  is  processed 
through  data  reduction  pipelines  to  measure  photometric  and  astrometric  prop¬ 
erties  (Stoughton  et  al.  2002;  Smith  et  al.  2002;  Pier  et  al.  2003;  Abazajian  et  al. 
2003)  and  to  identify  targets  for  spectroscopic  follow-up. 


where  /?  =  7/rcut  ^  1/24  in  our  calculations.  To  constmct  the 
flattened  (oblate)  DMHs,  a  nonhomologous  transformation  is 
applied  to  equation  (3)  to  achieve  the  desired  axis  ratio  qh  while 
preserving  the  central  density.  In  order  to  minimize  computa¬ 
tional  time  when  constmcting  flattened  DMHs  with  embedded 
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TABLE  3 
Parameter  Space 


Best  Fits 


Prograde  Retrograde 

l-b  Range  No.  of  -  - 

Parameter  (deg)  Precision  Values  prol  pro  2  pro  3  retl  ret2  ret 3 


Ms  (xlO8  Mq) .  [0.6,12.0]  ±3.0  3  6.3  6.3  6.3  12.0  12.0  12.0 

ra  (kpc) .  [17.5,80.5]  ±3.5  10  22.8  22.8  22.8  22.8  22.8  22.8 

e .  [0,  0.7]  ±0.05  7  0.1  0.1  0.1  0.5  0.5  0.5 

i  (deg) .  [5.0,  45.0],  [135.0,  175.0]  ±5.0  10  25.0  25.0  25.0  165.0  155.0  165 

(j)  (deg) .  [0.0,360.0]  ±15  12  300  150.0  120.0  180  210  210 

qh .  [0.0,360.0]  ±0.05  6  0.6  0.7  0.8  0.6  0.6  0.5 

Xbest .  6.3  ...  ...  1.92  1.96  1.97  2.34  2.39  2.39 


Note. — All  quantities  are  given  at  t  =  0  of  the  simulation;  ra  denotes  the  initial  distance  to  the  Galaxy  center;  e  =  (ra  —  rp)/(ra  ±  rp)  is  the  orbital  eccentricity;  ra 
and  rp  are  the  apo-  and  pericenter  distances,  respectively;  i  is  the  inclination  with  respect  to  the  disk  plane;  <f>  is  the  azimuthal  angle  in  spherical  coordinates;  and  qh  is 
the  axis  ratio  of  the  halo’s  density  profile.  Values  of  free  parameters  are  equally  distributed  with  the  ranges.  Our  notation  is  so  that  0°  <  i  <  90°  indicates  a  prograde 
motion,  whereas  90°  <  i  <  180°  a  retrograde  one. 


bulges  and  disks,  we  apply  a  highly  efficient  technique  using 
multipole  potential  expansions  to  tailor  the  local  velocity  ellipsoid 
to  the  required  morphology  (Body  et  al.  2001).  The  algorithm  to 
add  together  individual  components  in  a  single  galaxy  is  adapted 
from  Hemquist  (1993).  This  code  scales  linearly  with  particle 
number  and  hence  we  can  constmct  flattened  DMHs  with  ^1 06 
particles  in  a  short  computational  time. 

The  specific  parameters  are  chosen  to  reproduce  the  observed 
rotational  curve  of  the  Milky  Way.  The  only  free  parameter  of  the 
Galaxy  model  that  we  explore  in  our  calculations  is  the  halo  axis 
ratio  (qh). 

The  dwarf  satellite  model  follows  a  King  profile  (King 
1966)  with  central  potential  TGr2  =  4  and  concentration  c  = 
log10(G/rK)  =  0.84,  where  rt  and  rK  are  the  tidal  and  King  radii, 
respectively.  Since  we  are  assuming  that  the  tidal  debris  repre¬ 
sent  multiple  wraps  (and  thus,  multiple  perigalacticon  passages) 
of  the  tidal  stream  originating  from  one  disrupting  satellite,  we 
restrict  the  tidal  radius  of  the  progenitor  dwarf  galaxy  in  each 
model  with  the  following  criteria:  (1)  the  tidal  radius  must  be 
small  enough  to  prevent  galaxy  disruption  by  the  first  peri¬ 
galacticon  passage,  and  (2)  the  tidal  radius  must  be  large  enough 
to  induce  a  progressive  mass  loss  that  leads  to  the  formation  of  a 
tidal  stream.  We  have  selected  the  satellite  tidal  radius  to  match  the 
Jacobi  limit  (§  7.3  of  Binney  &  Tremaine  1987)  at  the  first  peri¬ 
galacticon,  which  is  a  simple  method  to  achieve  “slow”  mass  loss 
along  the  orbit. 

3.2.  Orbital  Parameters 

Given  the  Milky  Way  and  satellite  mass  profiles,  three  aspects 
determine  the  subsequent  orbit  evolution:  (1)  the  flattening  of  the 
Milky  Way  halo  qh;  (2)  the  mass  of  the  satellite,  determining  the 
degree  of  dynamical  friction;  and  (3)  the  initial  orbital  geometry. 
The  corresponding  free  parameters  are  outlined  in  Table  3.  We 
have  performed  151,200  simulations  by  means  of  our  semi- 
analytic  code  (see  below)  in  order  to  fit  the  observational  data. 
Each  orbit  is  repeated  with  12  different  azimuthal  angles  since,  a 
priori,  we  do  not  know  whether  dynamical  friction  affects  the  orbit 
within  one  orbital  period  (which  would  “break”  the  axisymmetry 
of  the  problem). 

Because  of  the  high  efficiency  of  the  semianalytic  code,  the 
resolution  achieved  in  the  parameter  space  is  constrained  by  the 
limited  number  of  available  observational  data  rather  than  by 
CPU  limitations.  More  densely  sampled  parameter  surveys  would 
not  provide  stronger  constraints,  owing  to  the  large  degeneracy 


seen  already  in  this  sparse  sampling  of  parameters  space  (see 
§4.1). 

3.3.  Semianalytic  Fit 

Exploring  a  large  number  of  initial  conditions  and  accounting 
for  the  time-dependent  dynamical  friction  and  satellite  disrup¬ 
tion  can  only  be  carried  out  in  a  reasonable  time  by  using  semi¬ 
analytic  algorithms.  We  have  used  that  proposed  by  Penarrubia 
(2003),  which  provides,  for  the  galaxy  density  profile  outlined  in 
§  3.1,  the  evolution  of  r,  where  r  is  the  satellite  center-of-mass 
position,  and  the  dwarf  satellite  mass  Ms.  This  code  has  been 
tested  against  A-body  calculations  for  a  large  spectrum  of  orbital 
parameters  and  satellite  masses  (Just  &  Penarrubia  2005),  as 
well  as  halo  flattenings  (Penarrubia  et  al.  2004),  showing  that 
I  ^analytic  —  **v-body|  <0.7  kpc  for  time  integrations  of  3  Gyr,  once 
the  Coulomb  logarithm  is  fit  to  A-body  simulations  (the  best-fit 
values  being  In  =  2.1  for  the  halo  and  In  =  0.5  for  the 
disk). 

In  order  to  compare  our  collection  of  orbits  to  the  observa¬ 
tional  data  we  assume  that  the  distribution  of  stream  stars  can  be 
reproduced  by  a  single  stream  within  a  given  number  of  wraps. 
Equivalently,  we  can  define  the  number  of  wraps  as  the  time  in¬ 
terval  that  a  point-mass  particle  needs  to  cover  the  phase  space 
defined  by  the  observational  data,  AT,  which  a  priori  is  an  un¬ 
known  quantity  that  we  have  estimated  from  our  A-body  simu¬ 
lations  to  lie  between  2Torb  <  AT  <  37orb,  where  Torb  =  2tt rlvc 
is  one  dynamical  period  at  r  ~  20  kpc  from  the  Galaxy  center; 
vc  =  220  km  s-1.  The  maximum  integration  time  of  semianalytic 
orbits  is  Tmax  =  5Torb  ~  3  Gyr. 

The  fitting  algorithm  that  we  use  is  the  following: 

1 .  Each  orbit  is  divided  into  n  overlapping  segments  of  ex¬ 
tension  AT  Following  the  above  discussion,  we  carry  out  the  fit 
for  stream  extensions  of  AT  =1.2  and  1.8  Gyr  (i.e.,  AT  ~  2 
and  3Torb).  The  number  of  segments  for  a  single  orbit  is  n  =  4 
and  3,  respectively.8 

2.  For  a  given  segment  we  calculate 

2  f  C,  analytic  A  obs^\  ^  f  Vj,  rad,  analytic  A  rad,obs^\ 

2  \  Or,  i  )  \  Ov  i  J 

8  This  selection  leads  to  fixed  overlapping  time  intervals  of  extension 
(nAT  —  Tmax)/(n  —  1),  i.e.,  0.6  Gyr  for  n  —  4  and  1.2  Gyr  for  n  =  3. 


132 


PENARRUBIA  ET  AL. 


Vol.  626 


where  vrad  is  the  heliocentric  radial  velocity  and  ar_  /  and  av p 
are,  respectively,  the  errors  in  the  distance  estimation  and  in 
the  radial  velocity  of  the  observational  point  i.  The  distance 
error  is  approximately  15%  of  the  heliocentric  distance,  whereas 
crv  ~  5  km  s_1.  The  subindexes  i  and  j  take  the  values  i  =  1, 
2,  .  .  .  ,  TVobs,  where  Aobs  is  the  number  of  observational  data,  j  = 
1,  2,  .  .  .  ^analytic?  and  ^analytic  is  the  number  of  semianalytic 
points  in  a  segment.  Subsequently,  for  each  observational  value 
i  we  look  for  the  semianalytic  point  that  leads  to  the  smallest 
value  of  Xij,  obtaining  a  set  of  Aobs  semianalytic  points  where 
the  x-values  find  a  minimum.  The  average  value  of  x  in  a  given 
segment  is  x  =  E  Xi^ohs. 

3.  The  calculation  is  repeated  for  the  rest  of  segments. 

4.  The  x-values  of  each  segment  are  sorted,  obtaining  the 
time  interval  of  the  orbit  where  the  best  fit  occurs  and  the  min¬ 
imum  x-  This  method  allows  us  to  fit  more  than  one  orbital 
period  to  observational  points  located  in  a  small  region  of  the 
space  that  were  possibly  stripped  in  different  orbital  wraps  (see 
§  4.2). 

5 .  The  process  repeats  for  the  whole  orbit  collection  (regarding 
that  each  orbit  corresponds  to  one  point  of  our  parameter  space). 

3.4.  N-body  Code 

The  semianalytic  code  reproduces  remarkably  well  the  motion 
of  a  dwarf  satellite  center  of  mass  and  its  mass  evolution.  How¬ 
ever,  as  Piatek  &  Pryor  (1995)  showed,  the  mass-loss  process 
itself  is  fairly  complicated  and  difficult  to  implement  in  semi¬ 
analytic  algorithms. 

In  order  to  describe  the  distribution  of  stripped  stars  in  the  Gal¬ 
axy  we  perform  A-body  simulations  from  the  best-fitting  orbits 
found  using  the  semianalytic  code.  We  carry  out  self-consistent 
A-body  calculations  with  SUPERBOX,  a  particle-mesh  code  (see 
Felhauer  et  al.  2000),  which  calculates  the  gravitational  potential 
in  three  boxes  centered  at  the  disk,  bulge,  halo,  and  dwarf  satellite, 
each  box  with  643  grid  cells.  We  refer  the  reader  to  Penarrubia 
(2003)  and  Penarrubia  et  al.  (2002)  for  a  detailed  description  of  the 
code  parameters.  Here  we  merely  comment  that,  after  fixing  the 
time  step  to  0.65  Myr,  we  obtain  a  conservation  of  total  energy  and 
total  angular  momentum  of  around  1  %. 

The  number  of  particles  of  each  subsystem  is  Nd  =  9.0  x  104 
(disk),A&  =  5.0xl03  (bulge),  A/*  =  1.2  x  106  (dark  matter  halo), 
and Ns  =  1.0x10s  (dwarf galaxy). 

As  in  the  case  of  semianalytic  calculations,  we  evolve  A-body 
satellites  approximately  3  Gyr.  Three  aims  lead  us  to  select 
this  integration  time:  (1)  to  reproduce  semianalytic  calculations 
thoroughly;  (2)  to  minimize  feedback  effects  from  the  host  gal¬ 
axy  affecting  real  dwarf  satellites  that  are  difficult  to  implement 
in  A-body  realizations  of  the  Milky  Way  (such  as  the  disk’s  spiral 
arms,  overdensity  regions,  and  warps;  those  effects  can  be  treated 
like  small  corrections  to  our  Galaxy  potential  and  thus  be  ne¬ 
glected  in  short-time  orbit  calculations);  and  (3)  to  approximate 
the  orbits  of  tidal  stream  stars  as  the  orbit  of  the  main  system. 
When  obtaining  the  progenitor’s  main  orbital  properties  from 
semianalytic  calculations,  one  assumes  implicitly  that  the  es¬ 
caping  particles  follow  the  main  system’s  orbit  during  a  given 
time  interval.  This  assumption  only  holds,  therefore,  for  a  limited 
number  of  orbital  periods. 

4.  RESULTS 

4.1.  Semianalytic  Solutions 

In  Figure  1  we  plot  the  results  of  our  semianalytic  fit  to  ob¬ 
servational  data  for  two  different  stream  extensions,  AT  ~  2  Torb 
and  3Torb,  and  for  prograde  and  retrograde  orbits,  as  the  sense  of 
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Fig.  1 . — Minimum  values  of  the  fitting  parameter  x  as  a  function  of  initial 
heliocentric  distance  {bottom  left  panel),  eccentricity  {bottom  right  panel),  halo 
axis  ratio  {top  right  panel)  and  orbital  inclination  {top  left  panel).  For  each 
parameter,  we  perform  the  fit  assuming  that  the  observed  tidal  stream  can  be 
reproduced  by  an  orbit  segment  of  extension  AT  ~  2Tor^  {solid  lines)  or  3  Tor^ 
{dotted  lines). 

rotation  has  not  yet  been  unambiguously  determined.9  This  fig¬ 
ure  shows  that  the  best-fitting  parameter  values  do  not  depend  on 
the  exact  value  of  AT  However,  larger  segments  do  alter  the  value 
of  the  minima  and  the  smoothness  of  the  curves  by  increasing  the 
ratio  A/Aobs  (i.e.,  in  a  given  segment,  more  semianalytic  points 
can  be  fit  to  the  same  number  of  observational  constraints),  im¬ 
proving  the  quality  of  the  fit  but  at  the  same  time  leading  to  less 
pronounced  minima,  which  hardens  the  selection  of  the  best-fit- 
ting  parameter  sample. 

In  the  bottom  right  panel  we  show  the  minimum  value  of  x 
as  a  function  of  the  orbit’s  initial  eccentricity.  The  minima  are 
located  at  e  ~  0. 1  for  prograde  orbits  and  e  =  0.5  for  retrograde 
ones.  As  we  see  in  the  bottom  left  panel,  the  initial  heliocentric 
distance  that  leads  to  best  fits  independent  of  the  orbital  sense  of 
motion  is  rs  ~  22.8  kpc. 

In  the  top  left  panel  we  show  the  dependence  of  Xmin  on  the 
inclination  angle  with  respect  to  the  disk  plane  (note  that  the 
orbital  inclination  of  retrograde  orbits  is  i  =  1 80  —  Ai  with  this 
notation).  We  find  that  the  orbital  inclination  of  retrograde  orbits 
cannot  be  accurately  determined.  These  orbits  reach  the  solar 
circle,  which  therefore  reduces  the  number  of  passages  by  the  Ga¬ 
lactic  anticenter  (where  most  of  the  observational  data  are  located). 
As  a  result,  retrograde  solutions  cannot  account  for  the  large 
vertical  dispersion  of  observational  points  (see  Fig.  3).  In  con¬ 
trast,  prograde  orbits  present  a  well-marked  minimum  at  i  ~  25° . 

The  absence  of  a  defined  vertical  structure  of  debris  also  leads 
to  a  degenerate  value  of  the  halo  flattening,  as  we  see  in  the  top 
right  panel.  For  prograde  orbits,  the  halo  axis  ratio  that  leads  to  the 
best  solution  is  for  qh  =  0.6,  although  similar  Xmin  values  have 
been  found  for  =0.7  and  0.8.  No  information  about  the  halo 
flattening  can  be  obtained  from  the  fit  if  orbits  are  retrograde. 

9  As  Martin  et  al.  (2004a)  pointed  out,  the  radial  velocity  curve  as  a  function  of  pro¬ 
jected  position  does  not  provide  by  itself  sufficient  information  to  determine  whether 
the  progenitor  of  the  Monoceros  stream  follows  a  prograde  or  a  retrograde  orbit. 
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Fig.  2. — Best-fitting  prograde  solutions  (see  Table  3).  Top  panel:  Heliocentric  distance  as  a  function  of  Galactic  longitude  (/).  Middle  panel:  Heliocentric  radial 
velocity.  Bottom  panel:  Projected  spatial  position  in  the  l-b  plane.  Open  triangles  and  filled  squares  represent  detections  of  stars.  Open  squares  denote  color-magnitude 
detections  (with  no  radial  velocity  measurements  available).  [See  the  electronic  edition  of  the  Journal  for  a  color  version  of  this  figure .] 


As  one  would  expect,  there  are  a  large  number  of  orbits  that 
lead  to  similar  values  of  %.  That  is  far  from  surprising  when  one 
takes  into  account  that  (1)  the  region  of  the  Galaxy  where  the 
Monoceros  stream  has  been  detected  is  relatively  small  (rse 
[12,  20]  kpc,  110°  </  <  240°,  in  Galactic  coordinates),  which 
represents  a  small  interval  of  the  orbit,  and  (2)  the  stream  presents 
a  large  dispersion  in  the  z-direction  (perpendicular  to  the  disk 
plane)  with  a  poorly  defined  structure  (see  Fig.  2). 

In  Table  3  we  summarize  the  fit  results  and  show  the  three 
best-fitting  retrograde  and  prograde  orbits.  The  values  of  xmin 
clearly  indicate  that  retrograde  orbits  lead  to  considerably  worse 

Xret  —  l-24xpro. 

Our  fit  technique  also  provides  a  coarse  estimate  of  the  pro¬ 
genitor’s  initial  mass,  which  we  estimate  to  be  around  6  x  108  M0 
if  it  moves  on  a  prograde  orbit  and  1  x  109  M0  if  the  orbit  is  ret¬ 
rograde.  We  remark  that  this  result  is  fairly  approximate  (see  pre¬ 
cision  estimates  in  Table  3),  since  the  orbital  decay  within  3  Gyr  is 
very  low  for  prograde  orbits  and  practically  negligible  for  retro¬ 
grade  ones  (note  that  the  effect  of  the  dwarf  satellite  mass  on  the 
semianalytic  orbit  calculation  occurs  through  dynamical  friction). 
The  poor  sensitivity  of  the  result  to  the  satellite  mass  justifies  our 
early  choice  of  a  coarse  grid  in  satellite  masses. 

In  Figure  2  we  compare  the  semianalytic  satellite  orbits  with 
the  observational  data  used  for  the  fit.  We  only  show  the  orbit  at 
the  time  interval  when  the  best  fit  to  the  observations  occurs  (see 


§  3.3).  The  best  prograde  orbits  stay  relatively  far  from  the  Gal¬ 
axy  center,  re  [20,  25]  kpc,  and  have  a  pericenter  at  /  ~  180°. 
Different  values  of  halo  flattening  and  initial  azimuthal  angle 
match  the  observed  heliocentric  distances,  radial  velocities,  and 
projected  positions  of  debris  with  a  similar  accuracy,  leading  to 
multiply  degenerated  solutions. 

In  Figure  3  we  illustrate  our  best-fitting  retrograde  orbits.  Those 
orbits  have  an  apogalacticon  at  /  ~  1 80°,  where  most  of  observa¬ 
tional  points  are  located.  They  reach  the  solar  circle,  moving 
within  r  G  [7.5,  22]  kpc.  Since  orbits  in  the  Galactic  potential  fol¬ 
low  rosettes,  the  retrograde  orbits  pass  the  Galactic  anticenter  con¬ 
siderably  less  frequently  than  the  prograde  ones. 

The  satellite  tidal  radius  was  selected  to  match  that  of  the 
Galaxy  at  the  pericenter;  therefore,  we  find  that  for  prograde  sat¬ 
ellites  (rK,  g)  =  (0.51,  1.79)  kpc  whereas  for  retrograde  ones 
(r k ,  rt)  =  (0 .08 ,  0 . 54)  kpc,  where  rK  and  r,  are  the  King  and  the 
tidal  radius,  respectively.  Note  that  our  retrograde  satellites  are 
approximately  3.3  times  smaller  than  the  progrades  in  order  to 
prevent  tidal  dismption  by  the  first  perigalacticon  passages. 

4.2.  N-body  Calculations:  Orbital  Sense  of  Motion 

In  this  section  we  analyze  the  resulting  kinematical  and  spatial 
distribution  of  debris  from  the  best-fitting  orbits  obtained  by  our 
TV-body  algorithm.  In  Figure  4  we  plot  the  projection  of  the  pro- 
and  retrograde  orbits  in  the  A-Tplane,  the  radial  velocity  curve, 
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Fig.  3. — Same  as  Fig.  2,  but  for  our  best-fitting  retrograde  solutions.  [ See  the  electronic  edition  of  the  Journal  for  a  color  version  of  this  figure .] 


and  proper  motions  (in  Galactic  longitudinal  and  latitudinal  com¬ 
ponents)  obtained  from  A-body  realizations  of  our  best-fit  models 
and  compare  them  to  the  observational  data  outlined  in  §  2.  We 
have  integrated  the  prograde  orbit  for  2.99  Gyr  in  order  to  repro¬ 
duce  the  observed  projected  distribution  of  debris  (see  §  4.3), 
whereas  the  retrograde  one  was  evolved  only  1 .9  Gyr,  until  it  was 
close  to  total  disruption.  Even  imposing  the  initial  satellite  tidal 
radius  to  match  the  Jacobi  limit  at  the  first  perigalacticon  distance 
(which  leads  to  a  retrograde  satellite  approximately  3.3  times 
smaller  than  the  prograde  one)  is  not  sufficient  to  assure  longer 
survival  times,  likely  because  of  the  enhanced  mass  loss  induced 
by  disk  and  bulge  tidal  shocks. 

4.2.1.  Debris  Kinematics 

As  noted  in  §  4. 1 ,  the  slope  of  the  radial  velocity  curve  nearby 
the  Anticenter  dvm&/dl,  l  =  180°  can  be  well  reproduced  either 
by  low  eccentricity,  prograde  orbits  with  perigalacticon  at  /  ~ 
1 80°  or  by  high  eccentricity,  retrograde  orbits  with  apogalacticon 
at  /  ~  180°.  Hence,  the  information  provided  by  the  radial  ve¬ 
locity  curve  is  not  sufficient  to  determine  the  rotational  sense  of 
motion  of  the  tidal  stream  progenitor.  We  must  remark  on  the 
notably  low  observational  velocity  dispersion,  which  is  better  re¬ 
produced  by  the  prograde  model,  as  one  would  expect  from 
semianalytical  results  (see  Fig.  2,  middle  panel). 

Proper  motions  in  the  latitudinal  direction  do  not  break  the 
degeneracy  between  the  two  different  orbital  senses  of  motion 
either,  since  both  models  predict  a  similar  range  of  pjj  values, 
consistent  with  the  observations  taking  into  account  the  large 


error  bars.  In  contrast,  proper-motion  constraints  in  the  Galactic 
plane  direction,  ph  settle  the  sense  of  rotation  to  be  prograde 
(Fig.  4,  right  panels,  bottom  plots).  The  prograde  model  predicts 
proper  motions  in  that  region  in  the  range  pt  G  (— 5 ,  0)  mas  yr-1, 
whereas  for  the  retrograde  model  we  find  a  much  larger  azimuthal 
angular  velocity  pt  e(5,  20)  mas  yr-1  due  to  the  larger  relative 
velocity  between  the  Sun  and  the  dwarf  galaxy.  In  order  to  con¬ 
clusively  differentiate  the  orbital  sense  of  motion,  we  use  proper 
motions  of  10  stars  in  the  stream  derived  by  Munn  et  al.  (2004). 
Comparing  these  measurements  with  the  theoretical  predictions 
we  find  that,  except  for  two  stars,  observed  values  of  pi  are  re¬ 
markably  well  reproduced  by  the  prograde  model  (in  magnitude 
as  well  as  in  sign). 

We  have  transformed  vr,  ph  and  p^  into  Cartesian  velocities  in 
the  Galaxy  frame  (eq.  [7])  in  order  to  determine  the  components 
of  angular  momentum  and  the  orbital  inclination  cos  i  =  — LJL , 
where  L 2  =  Zj  +  is  the  total  angular  momentum  per  unit 
mass.10  In  Figure  5  we  plot  the  distribution  of  A-body  particles 
(models  prol  and  retl)  and  observational  points  in  the  angular 
momentum  plane  LR-LZ  {top  panel).  Tidal  stream  particles  are 
located  in  well-defined  regions  determined  by  the  main  orbital 
inclination  {Lz  <  0,  Lz  >  0  for  prograde  and  retrograde  orbits, 
respectively)  and  by  the  main  eccentricity  (which  decreases  for 

10  We  note  that  the  observational  values  of  LR  are  positive  definite  quantities 
with  large  errors,  which  leads  to  a  biased  estimate  unless  some  statistical  cor¬ 
rection  is  applied.  We  have  corrected  them  by  using  the  technique  of  Wardle  & 
Kronenberg  (1974). 
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Fig.  4. — Prograde  (model  prol )  vs.  retrograde  (model  retl )  orbits.  Left  panels:  X-Y  Galaxy  plane  projections.  The  sun  is  placed  at  {X,  7,  Z)  =  (—8,  0,  0)  kpc. 
Right  panels:  Heliocentric  radial  velocity  curve  {top  plots;  given  in  km  s-1).  Proper  motions  in  the  latitudinal  {middle  plots )  and  in  the  longitudinal  {bottom  plots) 
components,  both  given  in  mas  yr-1.  [See  the  electronic  edition  of  the  Journal  for  a  color  version  of  this  figure .] 
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Fig.  5. — Top  panel:  Prograde  (model  prol )  vs.  retrograde  (model  retl)  values  of 
the  angular  momentum  components.  Dots  represent  observational  values  after  be¬ 
ing  corrected  (see  text).  Bottom  panel:  Orbital  inclination  as  a  function  of  Galactic 
longitude. 


decreasing  values  of  LR ).  In  the  bottom  panel  we  plot  the  orbital 
inclination  as  a  function  of  Galactic  longitude  (/  <  90°  prograde 
orbits,  i  >  90°  retrograde  orbits),  which  shows  that  the  Mono- 
ceros  stream  progenitor  likely  follows  a  prograde  orbit.  Unfor¬ 
tunately,  the  observational  errors  are  too  large  to  provide  an 
estimation  of  the  main  orbital  inclination  from  these  stars. 

It  is  interesting  to  note  that  simply  by  measuring  kinematical 
properties  of  stream  stars  (radial  velocities  plus  proper  motions) 
we  would  be  able  to  determine  not  only  the  sense  of  motion  of  the 
progenitor  system,  but  also  its  orbital  eccentricity  and  inclination 
if  those  measurements  were  accurate  enough. 

4.2.2.  Geometrical  Distribution  of  Debris 

As  Figure  5  shows,  the  kinematical  properties  of  tidal  stream 
particles  are  similar  to  those  of  the  progenitor  during  several  pe¬ 
riods.  As  a  result,  the  spatial  distribution  of  debris  is  an  approxi¬ 
mate  representation  of  the  dwarf  galaxy  orbit.  In  Figure  4  we  show 
the  projection  in  the  Galactic  X-Y plane  of  dwarf  galaxy  A-body 
particles  for  models  prol  and  retl ,  observing  the  following: 

1.  Stripped  particles  from  the  prograde  satellite  are  located 
preferentially  at  two  Galactocentric  distances,  r  ~  12  and  r  ~ 
22  kpc,  building  two  nearly  concentric  “rings”  in  space.  That 
peculiar  spacial  distribution  forms  because  of  ( a )  the  low  orbital 
eccentricity  and  (fe)  the  anisotropic  mass  loss.  As  Piatek  &  Pryor 
(1995)  found,  disrupting  dwarf  galaxies  become  elongated  ow¬ 
ing  to  the  action  of  tidal  fields.  Stripped  particles  escape  pref¬ 
erentially  parallel  to  the  major  axis,  which,  at  the  same  time,  is 
oriented  perpendicular  to  the  density  gradient  vector.  Because  of 
this  process,  particles  belonging  to  the  trailing  and  the  leading 
tails  move,  respectively,  outward  and  inward  with  respect  to  the 
main-body  orbit.  This  process  can  be  clearly  seen  in  the  top  left 
panel,  the  satellite  center  of  density  being  located  at  (A,  Y)  = 
(—15,  —15)  kpc  and  moving  clockwise. 

2.  The  retrograde  orbit  forms  a  rosette  with  apogalacticon  at 
22  kpc  and  perigalacticon  at  7.3  kpc,  also  providing  a  reasonably 


good  fit  to  the  Monoceros  star  detections  but  failing  to  reproduce 
the  gap  between  the  two  tails  observed  in  the  Tri/And  regions. 

4.3.  Position  of  the  Main  Body 

There  are  some  facts  that  point  to  the  possible  survival  of  the 
Monoceros  stream’s  progenitor.  Stars  escaping  from  a  tidally 
disrupting  system  follow  the  main  system  orbit  for  some  time 
(the  exact  evolution  of  their  orbits  depend  on  several  parame¬ 
ters),  building  up  what  has  been  defined  as  “tidal  tails.”  Stars 
that  were  lost  several  orbital  periods  ago  follow  orbits  in  the  host 
galaxy  potential  nearly  independent  of  their  parent  system.  These 
particles  would  spread  in  a  large  volume  of  space  (no  overdensity 
signature)  with  no  velocity  gradient.  Helmi  et  al.  (2003)  find  that 
only  careful  measurements  of  the  phase-space  structure  might 
indicate  whether  they  belonged  originally  to  a  disrupted  body. 
According  to  their  predictions,  the  spatial  location  of  debris  form¬ 
ing  a  typical  tidal  stream  structure,  as  well  as  their  well-defined 
radial  velocity  curve,  indicates  that  the  tidal  tails  that  we  observe 
are  possibly  young  and,  therefore,  either  (1)  the  main  body  has 
not  been  yet  completely  destroyed  or  (2)  the  dismption  occurred 
recently. 

As  shown  in  §  4.1,  the  range  of  distances  and  Galactic  latitudes 
where  parts  of  the  stream  are  detected  provides  robust  values  of 
the  orbital  inclination  and  eccentricity  of  the  progenitor.  However, 
the  exact  vertical  and  radial  distribution  of  debris  depends  on  the 
location  of  the  main  body  (or,  equivalently,  the  time  when  we  stop 
the  A-body  simulation). 

We  have  two  main  constraints  to  fix  the  progenitor  location: 
(1)  geometrical  distribution  of  debris  and  (2)  different  metal- 
licities  observed  in  different  regions. 

Geometrical  constraint. — Detections  obtained  so  far  show 
that  the  Monoceros  stream  forms  a  complex  vertical  struc¬ 
ture  from  the  disk  plane,  if  all  detections  belong  to  a  single  tidal 
stream.  In  this  paper  we  explore  the  scenario  in  which  such 
structure  arises  from  a  single  disrupted  satellite.  Different  studies 
(Majewski  et  al.  2004;  H.  J.  Newberg  et  al.  2005,  in  preparation) 
have  reported  the  presence  of  two  parts  of  the  stream  at  /  ~  110°, 
both  located  at  similar  projected  position  but  at  different  heliocen¬ 
tric  distances.  Observations  close  to  the  disk  plane  do  not  show  the 
presence  of  the  distant  tail  (fee  [—20°,  20°]),  likely  because  of 
dust  absorption,  up  to  b  >  20°,  where  it  is  again  detected.  We 
note  that  the  close  tail  has  been  detected  only  at  negative  latitudes 
(fee  [-20°,  -10°]). 

We  mainly  use  observational  points  at  1 10°  <  l  <  130°  to  de¬ 
termine  the  progenitor  position,  since  the  range  of  distances 
of  detections  at  /e  [140°,  240°]  is  considerably  smaller,  which 
therefore  provides  a  weaker  constraint.  Looking  at  this  longi¬ 
tudinal  range  we  observe  that  if  the  main  body  is  placed  at  /  > 
100°  (integration  time  t  <  2.66  Gyr),  we  find  no  distant  tail  at 
—25°  <  fe  <  —20°  (only  leading-tail  particles  moving  in  the  clos¬ 
est  “ring”  can  be  found  in  this  region),  whereas  if  we  place  the 
main  body  at  /  <  200°  (integration  time  t  >  3.04  Gyr),  we  find  a 
luminous  close  tail  at  fe  >  20°,  not  present  in  observations.  There¬ 
fore,  we  obtain  comparable  fits  to  observations  if  the  main  body  is 
located  in  the  range  100°  >l>  200°,  which  corresponds  to  in¬ 
tegration  times  from  2.66  to  3.04  Gyr,  respectively  (note  that  pro¬ 
grade  orbits  move  with  dll dt  <  0). 

Metallicity  constraint. — Detections  of  the  distant  tail  at  100°  < 
/  <  150°  show  metallicities  considerably  lower  than  those  of  the 
close  tail  in  the  same  longitudinal  range  (Rocha-Pinto  et  al.  2003). 
Although  we  have  no  information  on  the  initial  metallicity  gra¬ 
dient  of  the  progenitor,  we  assume  that  (1)  stars  with  originally 
high  binding  energy  are  those  with  the  highest  metallicity,  as  star 
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Fig.  6. — Distributions  of  initial  binding  energies.  We  compare  the  initial 
distribution  {solid  lines)  with  that  at  t  =  2.66  Gyr  {top panel)  and  t  =  2.99  Gyr 
{bottom panel)  for  particles  in  the  range  100°  <  /  <  150°.  For  each  integration 
time  we  distinguish  between  particles  in  the  “inner  ring”  (r  <  18  kpc  from  the 
Galaxy  center;  dotted  lines)  and  in  the  “outer  ring”(r  >18  kpc;  dashed  lines). 
At  t  =  2.66  and  t  =  2.99  Gyr  the  main  body  is  located  at  /mb  =  95°  and  /mb  = 
245°,  respectively.  We  assume  that  particles  with  initial  high  binding  energy 
\E!Emm\  ~  1  would  present  higher  metallicities  than  those  with  low  binding 
energy  at  t  =  0,  \E/Emm\  ~  0. 

formation  and  metallicity  enrichment  happen  more  intensely  in 
the  center  of  dwarf  galaxies  (Pagel  &  Edmunds  1981;  Harbeck 
et  al.  2001)  and  (2)  those  processes  generating  or  enhancing  any 
metallicity  gradient  have  stopped,  or  are  very  slow,  in  the  last 
3  Gyr.11 

Under  those  assumptions,  the  actual  distribution  of  stream 
metallicities  is  related  to  the  age  of  the  tidal  tails,  since  low- 
metallicity  stars  with  low  binding  energy  are  stripped  out  by  the 
action  of  tidal  forces  more  rapidly  than  stars  with  initially  high 
binding  energy  and  high  metallicity,  which  will  be  preferentially 
located  in  recently  stripped  tails. 

Matching  the  metallicity  observed  in  different  regions  of  the 
stream  reduces  the  range  of  possible  progenitor  locations  ob¬ 
tained  from  geometrical  and  kinematical  constraints.  We  have 
compared  relative  metallicities  between  the  distant  and  the  close 
tails  within  100°  <  /  <  150°.  In  Figure  6  we  plot  the  distribution 
of  stream  particles  as  a  function  of  their  initial  binding  energy  for 
two  different  positions  of  the  main  body  (equivalently,  two  final 
integration  times).  Solid  lines  show  the  initial  distribution  (once 
the  satellite  is  numerically  relaxed),  and  dotted  and  dashed  lines 
show  the  distribution  of  particles  with  r  <  1 8  kpc  (“inner  ring”) 
and  r  >  18  kpc  (“outer  ring”),  respectively. 

The  top  panel  of  Figure  6  shows  the  result  of  placing  the  sat¬ 
ellite  at  /  =  95°.  We  can  see  that  the  number  of  particles  with  low 
initial  binding  energy  (which  one  would  expect  to  be  metal-poor) 
is  larger  in  the  close  tail,  contradicting  observations.  If  we  inte- 


1 1  This  assumption  is  made  in  order  to  compare  consistently  the  initial  and  the 
final  distribution  of  binding  energies.  New  episodes  of  star  formation  during  the 
tidal  disruption  of  the  dwarf  galaxy  would  likely  occur  in  the  central  regions  of 
the  system  (thus,  increasing  the  metallicity  gradient),  which  does  not  alter  the 
conclusions  of  our  comparison. 


grate  a  longer  time,  so  that  the  main  body  locates  at  /  =  245°,  the 
number  of  particles  in  the  distant  tail  with  initial  low  binding  en¬ 
ergies  is  clearly  larger  than  in  the  close  tail.  In  this  case,  the  outer 
tail  has  in  average  lower  metallicity  than  the  close  one. 

It  is  interesting  to  note  that  the  ideal  mass-loss  process,  in 
which  mass  shells  are  progressively  removed,  is  fairly  approxi¬ 
mate.  Shell  crossing  occurs  during  the  evolution  of  the  satellite, 
likely  because  of  the  action  of  Galactic  tidal  forces  and  shocks. 
As  a  result,  we  observe  the  presence  of  particles  with  initially 
high  binding  energy  in  both  stream  tails.  This  fact  might  explain 
observations  of  high-  and  low-metallicity  stars  in  same  fields 
(e.g.,  Rocha-Pinto  et  al.  2003).  Although  some  mixing  occurs, 
Figure  6  shows  that  the  number  of  particles  with  initial  high  bind¬ 
ing  energy  that  remain  in  the  main  system  is  larger  than  in  the  tidal 
streams. 

In  Figure  7  we  plot  the  heliocentric  distance  {top  panel ) 
and  Galactic  latitude  {middle  panel)  of  debris  as  a  function  of 
Galactic  longitude.  The  bottom  panels  show  the  distribution  of 
particles  around  the  main  system  position.  The  combination 
of  geometrical  and  metallicity  constraints  fixes  the  main  body  at 
/  ~  245°  zb  3°,  b  ~  — 18°  zt  2°,  corresponding  to  an  integration 
time  of*  =  2.99  Gyr.  The  heliocentric  distance  is  approximately 
15.2  dz  0.8  kpc. 

We  note  that  this  result  is  not  definitive.  Unfortunately,  the 
available  data  are  not  yet  sufficient  to  provide  strong  constraints 
on  the  main  system  location.  The  progenitor  position  presented 
here  is  fairly  approximate  and  is  based  on  the  assumption  of  an 
initial  metallicity  gradient  of  the  progenitor.  We  expect  that  fu¬ 
ture  detections,  mapping  larger  areas  of  the  sky,  will  provide 
stronger  geometrical  constraints  and  reduce  the  number  of  pos¬ 
sible  scenarios. 

5.  DISCUSSION 

In  this  section  we  discuss  the  possible  identification  of  the 
Monoceros  stream  progenitor  in  CMa  (Martin  et  al.  2004a),  by 
comparing  its  orbital  properties  to  those  of  our  model  prol.  We 
also  analyze  the  orbits  of  three  globular  clusters  with  measured 
proper  motions  thought  to  be  associated  with  the  CMa  dwarf. 
Finally,  we  briefly  summarize  and  discuss  alternative  explana¬ 
tions  for  the  Monoceros  tidal  stream  that  can  be  found  in  the 
literature. 

5.1.  The  Triangulus /Andromeda  Stream 

In  Figure  8  we  show  the  location  of  the  recent  detected  Tri/And 
tidal  streams  (Rocha-Pinto  et  al.  filled  squares )  against  the 

debris  distribution  of  model  prol  {gray).  We  also  show  those  de¬ 
tections  used  to  constrain  the  orbital  parameters.  The  location  and 
radial  velocity  of  CMa  overdensity  region  is  shown  by  a  large 
open  circle.  The  prograde  model  suggests  the  recently  discovered 
stream  in  Triangulus/Andromeda  as  natural  part  of  the  Monoceros 
stream,  both  fitting  accurately  to  the  modeled  kinematics  and 
spatial  distribution  of  debris.  The  observed  Tri/And  streams  (/  G 
[120°,  150°])  appear  as  a  connection  between  the  Monoceros 
stream  (/  G  [160° ,  240°])  and  the  new  SDSS  and  Ibata  detections 
at  /  G  [  1 1 0° ,  1 40°  ] ,  justifying  the  inclusion  of  these  data  in  our  fit, 
as  we  discuss  in  §  2. 

5.2.  The  Canis  Major  Dwarf  as  the  Progenitor 
of  the  Monoceros  Tidal  Stream 

Martin  et  al.  (2004a)  identified  an  elliptical  overdensity  of 
M-giant  stars  in  the  CMa  region  with  properties  suggestive  of  a 
disrupted  dwarf  galaxy:  (1)  standard  models  of  the  Milky  way 
cannot  account  for  the  large  number  of  red  giants  in  that  region  of 
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Fig.  7. — Model  prol  integrated  for  2.99  Gyr.  Heliocentric  distance  {top  panel )  and  Galactic  latitude  {middle  panel)  of  debris  as  a  function  of  Galactic  longitude.  We 
use  notation  of  Fig.  4  to  distinguish  different  observational  sources.  Gray  and  black  dots  denote,  respectively,  particles  at  r  >  18  kpc  (“outer  ring”)  and  r  <  18  kpc 
(“inner  ring”)  from  the  Galaxy  center.  The  bottom  panels  show  the  distribution  of  debris  around  the  main  system  remnants.  From  left  to  right  we  plot  the  heliocentric 
distance  and  the  Galactic  longitude  and  latitude  distributions,  which  are  normalized  to  the  total  number  of  dwarf  galaxy  stars  {Ns  =  105). 
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our  Galaxy,  and  (2)  the  low  dispersion  of  the  radial  velocity  dis¬ 
tribution  (approximately  20-25  km  s-1)  is  unexpected  for  a  disk 
population.  From  the  number  of  M-giant  stars,  these  authors  es¬ 
timate  a  mass  of  108— 109  M0,  similar  to  the  Sagittarius  dwarf 
galaxy.  The  location  in  Galactic  coordinates  is  220°  <  /  <  260°, 
— 15°  <  b  <  —7°,  with  no  data  for  b  >  —7°  because  of  dust  ab¬ 
sorption.  In  a  second  paper,  Martin  et  al.  (2004b)  find  that  the 
heliocentric  distance  of  this  systems  is  rs  ~  7.2  zb  2  kpc,  with  the 
maximum  surface  density  being  located  at  /  =  240°  and  b  = 
— 8?8 .  They  also  measure  a  the  radial  velocity  of  vr  =  1 09  km  s“ 1 . 
Adopting  the  selection  criteria  of  Bellazzini  et  al.  (2004)  and 
Martin  et  al.  (2004a)  to  select  M-giant  stars  from  the  CMa  field, 
Momany  et  al.  (2004)  measure  the  following  proper  motions 
(w,  M*)  =  ( -  3-5  ±  2,  -o.l  ±  2)  mas  yr_1. 

Martin  et  al.  (2004a)  concluded  that  CMa  is  a  satellite  galaxy 
undergoing  tidal  disruption,  and  because  of  its  apparent  similar 
position,  they  suggested  that  it  probably  is  the  remnant  of  the 
Monoceros  stream  progenitor.  Here  we  compare  the  orbital  prop¬ 
erties  of  CMa  to  our  best  fit  of  the  Monoceros  stream  progenitor 
in  order  to  analyze  a  possible  common  origin  of  both  systems. 


We  must  remark  that  our  comparison  is  still  fairly  preliminary 
since  it  goes  beyond  the  original  scope  of  this  paper. 

5.2.1.  Orbit  Calculation 

We  have  integrated  the  orbit  of  CMa  back  in  time  in  order  to 
compare  its  kinematical  properties  with  those  of  our  model.  The 
velocity  of  CMa  in  Cartesian  coordinates  with  origin  in  the 
Galactic  center  was  obtained  using  the  following  expression: 

r  =  r0  +  i;rad(cos  b  cos  /,  cos  b  sin  /,  sin  b) 

"T  — cos  b  sin  /,  cos  b  cos/,  0) 

-f  r^/^— sin  Z?  cos /,  —  sin b  sin/,  cosZ?),  (6) 

where  r0  =  (10.0,  225.2,  7.2)  km  s_1  (Binney  &  Merrifield 
1998)  is  the  solar  velocity,  rs  the  heliocentric  distance,  and  k,  ~ 
4.74  a  conversion  factor  from  (kpc  mas  yr_1)  to  (km  s-1).  The 
velocity  vector  is  (—145.4,  193.2,  —4.2)  km  s-1.  The  result¬ 
ing  CMa  orbit  has  low  inclination  (/  =  4°)  and  is  nearly  circular 
e  ~  0.16,  as  argued  by  Momany  et  al.  (2004).  Figure  9  shows 
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Fig.  8. — Comparison  of  our  model  for  the  Monoceros  stream  progenitor  against  CMa  dwarf  properties  (open  circles ).  Top  panel:  Heliocentric  distances  of  stream 
particles  from  model  prol  (gray)  against  observations.  Middle  panel:  Projection  in  Galactic  coordinates.  Bottom  panel:  Heliocentric  radial  velocities. 


Fig.  9. — Orbital  inclination  of  CMa  as  a  function  of  pi,  p^,  fixing  the  radial 
velocity  to  109  km  s-1  and  the  heliocentric  distance  to  rs  =  7.2  kpc. 


that  the  value  of  i  is  strongly  sensitive  to  the  errors  in  proper  mo¬ 
tions.  Taking  into  account  the  errors  estimated  by  these  authors, 
we  obtain  that,  within  1  a  error,  the  orbital  inclination  of  CMa 
lies  in  the  range  /  =  [0°,  1 8°  ]  if  we  fix  the  radial  velocity  to  vr  = 
109  km  s-1  and  the  heliocentric  distance  to  7.2  kpc.  Com¬ 
paring  this  result  with  the  best-fitting  orbits  that  we  present  in  this 
work,  one  can  see  the  following: 

1 .  The  orbital  eccentricity  is  fairly  similar  to  the  one  we  find 
(model prol  predicts  e  =  0.10  =b  0.05). 

2.  The  orbital  inclination  of  CMa  is  too  low  to  account  for  the 
large  vertical  dispersion  of  stream  stars.  This  result  is,  however, 
not  conclusive,  since  a  orbit  similar  to  that  of  model  prol  lies 
within  1.5  a  error.  More  thorough  measurements  of  proper  mo¬ 
tions  are  necessary  to  clarify  the  orbit  of  the  CMa  dwarf. 

3.  The  geometrical  and  metallicity  constrains  that  we  can 
impose  with  the  available  observational  data  fix  the  final  position 
of  our  best-fitting  model  at  /  ~  245°,  b  ~  —  1 5°  and  a  distance  to 
the  Sun  of  rs  ~  1 5  kpc.  CMa  appears  in  a  direction  similar  to  that 
of  the  main  system  of  model  prol ,  but  7  kpc  closer.  If  CMa 
proves  to  be  the  Monoceros  stream  progenitor,  a  way  to  reconcile 
the  close  distance  of  the  CMa  with  the  distant  stream  detections 
might  be  found  by  increasing  the  progenitor  mass.  A  simple 
estimation  of  the  decay  rate  induced  by  dynamical  friction  on 
nearly  circular  orbits  is  A r  oc  MsrAt  (eqs.  [7]— [25]  of  Binney  & 
Tremaine  1987).  For  model  prol ,  Ar  ~  3  kpc  after  At  =  3  Gyr. 
Imposing  A r  ~  9  kpc  in  3  Gyr  leads  to  an  initial  mass  of  the 
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progenitor  of  approximately  M'/Ms  =  Ar'/Ar  =  9/3  =  3,  i.e., 
M'  ~  1.8  x  109  M®.  However,  it  is  unclear  whether  the  tidal 
debris  of  such  a  massive  satellite  would  also  reproduce  the  ob¬ 
served  distribution.  In  order  to  check  this  assumption,  new  in¬ 
vestigations  adopting  the  CMa  dwarf  as  the  Monoceros  stream 
progenitor  appear  necessary,  which  goes  beyond  the  scope  of 
this  paper. 

5.2.2.  Three-dimensional  Velocity  Distribution 
in  the  CMa  Central  Region 

We  have  also  analyzed  the  distribution  of  radial  velocities  and 
proper  motions  in  different  regions  of  the  stream  in  order  to 
contrast  them  to  the  recently  available  observational  values. 

As  we  can  observe  in  Figure  8,  our  model  predicts  a  hot  tail 
overlapping  the  progenitor  remnants.  This  tail  was  stripped  out 
approximately  0.6  Gyr  ago  and  belong  to  the  trailing  tail  (i.e.,  to 
the  “inner  ring”).  In  the  top  left  panel  of  Figure  10  we  plot  the 
radial  velocity  distribution  (< dotted  line)  of  50  particles  located  in 
the  neighborhood  of  the  main  system’s  projected  position  and 
compare  it  to  that  of  the  center  of  our  satellite  model  {solid  line ; 
note  that  each  distribution  is  separately  normalized  to  the  number 
of  particles  in  the  sample).  This  figure  shows  than  the  velocity 
distribution  of  both  tails  can  be  clearly  differentiated,  with  peaks 
at  vr  ~  135  km  s-1  (cold  tail)  and  vr  ~  75  km  s-1  (hot  tail).  Al¬ 
though  some  particles  of  the  underlying  tail  are  included  in  the 
velocity  distribution  of  the  satellite  center,  one  cannot  identify  a 
second  peak  because  the  surface  density  of  the  hot  tail  is  much 
lower  that  of  the  satellite  center.  In  the  bottom  left  panel  of  this 
figure,  we  repeat  the  calculations  in  a  crossing  tail  region,  where 
the  surface  densities  of  both  tails  are  comparable.  In  this  case, 
one  observes  clearly  two  distributions  of  radial  velocities  (those 
from  the  distant  and  the  close  tail)  in  the  same  particle  sample. 

The  observed  low-dispersion  bimodality  cannot  be  reproduced 
by  standard  models  of  the  Milky  Way  and  is  likely  a  proof  of  a 
stream  detection,  as  we  show  here. 

Intriguingly,  Martin  et  al.  (2004b)  report  a  bimodal  velocity 
distribution  in  CMa  stars  located  at  the  center  of  the  CMa  dwarf. 
Measuring  the  radial  velocity  of  27  M-giant  stars  in  1  °  radius 
around  CMa  position,  they  observe  two  peaks  in  the  radial 
velocity  at  vr  ~  63  km  s_1  (10  stars)  and  vr  ~  109  km  s_1 
(17  stars),  with  very  low  dispersion  (around  5  and  11  km  s-1, 
respectively).  These  values  are  around  20%  larger  than  those  of 
model  prol.  Extrapolating  the  resulting  distribution  of  our  model 
to  their  detections  suggests  that  in  the  CMa  region  two  stream 
tails  are  overlapping.  Moreover,  the  fact  Martin  et  al.  (2004b) 
measure  similar  number  of  stars  in  both  peaks  indicates  that  their 
surface  density  must  be  similar  (if  we  assume  that  the  number  of 
M  giants  in  a  given  stream  is  proportional  to  its  surface  density). 
This,  however,  contradicts  the  presence  of  a  dwarf  galaxy  in 
Canis  Major  since  the  surface  density  at  the  center  of  these 
objects  is  several  orders  of  magnitude  larger  than  in  the  tidal 
streams. 

That  paradox  might  be  solved  if  their  small  sample  of  M-giant 
stars  were  strongly  contaminated  with  stars  that  belong  to  the  hot 
tail  surrounding  CMa.  In  that  case,  the  proper  motions  provided 
by  these  authors  might  be  also  affected.  In  right  column  of  Fig¬ 
ure  10  we  show  the  distribution  of  proper  motions  in  the  regions 
indicated  above.  As  in  the  distribution  of  radial  velocities,  dif¬ 
ferent  stream  parts  lead  to  well-differentiated  curves.  Looking 
at  the  region  where  the  dwarf  remnants  are  located,  we  can  ob¬ 
serve  that  the  main  system  shows  a  fairly  narrow  distribution, 
with  maxima  at  ((/Uz),  (pb))  —  (—1.7,  —0.6)  mas  yr-1,  whereas 
the  distribution  of  the  hot,  overlapping  stream  tail  is  centered  at 
((///),  {pb))  ~  (—  3.4,  3.4)  mas  yr-1.  Therefore,  including  hot 
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Fig.  10. — Left  panels:  Distribution  of  radial  velocities  in  the  CMa  region 
from  model  prol  (top)  and  in  a  region  of  tail  overlap  (bottom).  The  velocity 
distribution  of  the  satellite  center  includes  particles  within  243°  <  /  <  247°  and 
—  17  <  b  <  —13°.  The  number  of  particles  in  each  distribution  was  8400  in  the 
dwarf  center  and  50  in  the  surrounding  tail.  Dotted  and  dashed  lines  in  the  bottom 
panel  show  the  distributions  from  particles  in  the  inner  and  outer  “rings,”  re¬ 
spectively.  Right  panels:  Distribution  of  proper  motions  in  the  regions  indicated 
above. 


tail  stars  in  the  dwarf  remnant  sample  (in  a  significant  propor¬ 
tion)  would  lead  to  a  smaller  (pi)  and  a  larger  (pb),  which  results 
in  a  lower  orbital  inclination  (see  Fig.  9). 

5.3.  Possible  Associated  Stellar  Clusters 

Frinchaboy  et  al.  (2004)  have  collected  a  set  of  15  globular 
and  open  clusters  that  show  a  trend  in  their  radial  velocity  curve, 
as  well  as  in  spatial  location,  which  may  indicate  that  those  sys¬ 
tems  were  stripped  from  a  satellite  in  a  disruption  process. 

Martin  et  al.  (2004a)  claimed  that  four  of  those  clusters  (NGC 
1851,  1 904,  2298,  and  2808)  belong  to  a  globular  cluster  system 
associated  with  the  CMa  dwarf,  arguing  for  the  possible  detec¬ 
tion  of  the  Monoceros  stream’s  progenitor  in  Canis  Major. 

We  have  earlier  shown  that  radial  velocities  and  positions  do 
not  provide  sufficient  information  to  distinguish  between  pro- 
and  retrograde  orbital  motions.  In  order  to  break  that  degeneracy, 
accurate  measurements  of  proper  motions  are  needed.  In  this 
section  we  test  the  possible  association  of  those  clusters  with 
measured  proper  motions  (NGC  1851,  1904,  and  2298  from 
Dinescu  et  al.  1999)  with  the  CMa  dwarf.  In  addition,  we  also  dis¬ 
cuss  the  possible  common  origin  of  six  “possibly  associated 
clusters”  suggested  by  Frinchaboy  et  al.  (2004)  by  contrasting  lo¬ 
cations  and  radial  velocities  with  the  predictions  of  our  best-fitting 
model.  Unfortunately,  no  proper  motions  are  available  for  those 
systems,  so  the  results  are  not  conclusive. 

In  Figure  11  we  plot  the  X-Y projection  of  model  prol  {top 
panel )  and  clusters  listed  to  the  right  of  the  panel.  In  the  bottom 
panel,  the  top  plot  shows  the  radial  velocity,  and  the  middle  and 
bottom  plots  show  the  proper  motions  in  the  longitudinal  (/i/)  and 
latitudinal  (fib)  directions,  respectively.  This  figure  shows  the  fol¬ 
lowing  with  respect  to  CMa  and  other  clusters: 
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Fig.  1 1 . — Top  panel:  X-7projection  of  the  model  prol  debris  against  some  globular  and  open  clusters  possibly  associated  with  the  Monoceros  stream.  Bottom  panel, 
top  plot:  Radial  velocity  curve.  Bottom  panel,  middle  plot:  Proper  motions  in  the  latitudinal  component.  Bottom  panel,  bottom  plot:  Proper  motions  in  the  longitudinal 
component.  The  open  circle  denotes  the  position  and  radial  velocity  of  the  CMa  dwarf  (Martin  et  al.  2004b). 


CMa  clusters. — None  of  these  globular  clusters  present  ki- 
nematical  properties  consistent  with  those  of  the  CMa  dwarf. 
First,  NGC  1851  and  NGC  1904  have  radial  velocities  of  ap¬ 
proximately  320  and  208  km  s_1,  whereas  that  of  the  CMa  dwarf 
is  approximately  109  km  s-1  (Martin  et  al.  2004b).  After  inte¬ 
grating  the  orbit,  we  find  that  both  globular  clusters  move  on 
orbits  with  much  higher  eccentricities,  e  =  0.7  and  e  =  0.65, 
respectively  (in  agreement  with  Dinescu  et  al.  1999,  who  use  a 
similar  Galaxy  potential).  On  the  other  hand,  the  orbital  eccen¬ 
tricity  of  escaping  particles  remains  fairly  similar  to  that  of  the 
main  body  (see  Fig.  5),  leading  to  a  well-defined  radial  veloc¬ 
ity  curve.  If  one  assumes  that  this  holds  for  stripped  globular 
clusters,  we  come  to  the  conclusion  that  NGC  1851  and  1904 
are  unlikely  associated  with  CMa.  In  addition,  the  orbital  in¬ 
clinations  of  those  systems  are  considerably  higher  than  that  of 
CMa. 


The  third  globular  cluster  with  measured  proper  motions, 
NGC  2298,  follows  a  retrograde,  highly  eccentric  orbit  ( e  ~  0.78, 
Dinescu  et  al.  1999),  arguing  against  any  association  with  CMa 
and  the  other  clusters. 

Unfortunately,  there  are  no  proper  motions  available  for  NGC 
2208,  and  therefore  there  is  insufficient  information  to  allow  a 
determination  of  its  possible  association  with  CMa. 

Other  clusters. — In  the  middle  plot  of  the  bottom  panel  of 
Figure  11  we  see  that  these  clusters  present  a  trend  in  radial 
velocities  shown  by  Frinchaboy  et  al.  (2004)  consistent  with  our 
model.  Their  projected  positions  {top panel )  also  appear  to  form 
a  “stream”  in  space.  Searching  through  the  Galactic  globular 
cluster  sample  provides  some  additional,  plausible  candidates  to 
be  associated  with  the  Monoceros  stream:  Rup  1 06,  IC  4499,  and 
NGC  6779,  from  a  comparison  of  their  projected  position,  dis¬ 
tance,  and  radial  velocity  with  the  predictions  of  our  model. 
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Unfortunately,  whether  or  not  these  clusters  are  associated  with 
the  Monoceros  stream  cannot  be  asserted  in  the  absence  of  proper 
motions. 

5.4.  Thick-Disk  Stars  or  Tidal  Stream  Debris? 

The  main  selection  criteria  used  to  discriminate  between  Milky 
Way  stars  and  stars  of  an  “external”  origin  were  that  (1)  stream 
stars  are  observed  to  be  in  overdense  regions  not  predicted  by 
our  standard  Milky  Way  model  and  (2)  the  turn-off  stars  appear 
blue,  old,  and  metal-poor,  characteristics  similar  to  those  of  the 
Sagittarius  dwarf  galaxy  stars  (Newberg  et  al.  2002).  Those  con¬ 
ditions  by  themselves  do  not  exclude  other  explanations  for  the 
observations. 

Here  would  like  to  comment  briefly  on  two  possible  origins 
of  the  observed  tidal  streams:  (1)  thick  disk  and  (2)  a  satellite 
disruption. 

In  §  4  we  have  shown  that  the  stream  detections  can  be  re¬ 
produced  by  a  disrupting  satellite  on  a  prograde,  low-inclination, 
nearly  circular  orbit.  The  kinematical  and  as  spatial  distribution 
of  debris  resembles  what  one  would  expect  for  thick-disk  stars 
(see  Fig.  4)  except  for  one  point:  it  is  difficult  to  reconcile  the 
large  vertical  dispersion  (Az  =  2  sin  i  x  24.5  kpc  ~  20  kpc)  and 
the  large  distances  (~20  kpc  from  the  Galaxy  center)  of  proposed 
stream  stars  with  a  disklike  distribution.  To  clarify  this  point,  we 
have  carried  out  a  simple  experiment:  we  have  calculated  the 
probability  of  finding  thick-disk  stars  at  the  Monoceros  stream 
location.  The  number  of  stars  between  m  and  m  +  dm ,  where  m  is 
apparent  magnitude,  can  be  written  as 


Fig.  12. — Isocontours  of  the  thick-disk  number  of  particles  against  stream  de¬ 
tections.  The  Sun  is  located  at  R  —  z  =  0.  The  strong  line  represents  the  maximum 
likelihood  directions  of  finding  thick-disk  stars  toward  the  Galactic  anticenter  (/  = 
1 80°).  Lines  show  the  R,  z  values  where  the  number  of  thick-disk  stars  is  10%,  50%, 
and  90%  of  the  total  number  if  looking  at  the  anticenter  (solid  lines )  and  70°  away 
(dotted  lines). 


N(m)dm  =  N(rs)^-drs  =  N(rs)  -  -  -  — ,  (7) 

drs  In  10  rs 

where  rs  is  the  distance  of  a  star  to  the  Sun.  We  have  used  the  re¬ 
lationship  m  —  M  =  5  log10(r5)  to  convert  distances  into  apparent 
magnitudes.  The  probability  of  finding  a  star  between  rs  and  rs  + 
drs,  in  the  range  of  solid  angles  U,  U  +  dkl  is 


N(rs,  U)  drs  dkl  =  p{rs ,  kl)r^  drs  dfl,  (8) 

where  p{rs ,  U)  is  the  thick-disk  density  distribution  in  the  solar 
frame.  Combining  both  equations,  one  has  that  the  probability 
function  is 


P(rs,ft)  =Ap(rs,  Q)rs,  (9) 


where  A  is  some  normalization  constant.  The  fraction  of  thick- 
disk  stars  in  a  given  solid  angle  is  therefore 


N(rs,  SI) 


fgS  p(x,  fl)  x  dx 
f0°°  p(x,  U)  x  dx 


(10) 


As  we  can  see,  P  is  zero  for  rs  =  0  and  rs  — >  oo,  so  that,  for  a 
given  solid  angle,  there  is  a  distance  where  the  probability  finds  a 
maximum.  In  Figure  12  we  plot  the  maximum  likelihood  iso¬ 
contour  of  observing  thick-disk  stars  {thick  solid  line)  in  the  Ga¬ 
lactic  anticenter  direction  (eq.  [9]).  The  Sun  is  placed  at  R  =  z  = 
0.  Solid  lines  show  the  isocontours  of  the  thick-disk  fraction 
of  particles  as  a  function  of  their  position  in  cylindrical  coordi¬ 
nates  when  looking  in  the  Galactic  anticenter  direction  (eq.  [10]). 
Dotted  lines  show  the  isocontours  if  looking  70°  away  (i.e.,  /  = 
1 10°  or  /  =  250°)  in  order  to  take  into  account  the  range  of  Ga¬ 
lactic  longitudes  where  the  stream  has  been  observed.  In  both 
cases  we  plot  the  10%,  50%,  and  90%  contours.  As  we  can  see, 
the  Galactic  anticenter  is  the  direction  for  which  the  number  of 


distant  thick-disk  stars  is  maximum.  For  this  plot,  we  have  used 
the  thick-disk  model  proposed  by  Chen  et  al.  (2001),  which  fol¬ 
lows  exponential  profiles  in  the  planar  and  vertical  directions 
with  scale-lengths  of  3.5  and  0.75  kpc,  respectively. 

This  figure  shows  that  stream  detections  presented  in  this  work 
are  unlikely  related  to  the  thick-disk  population  because  (1)  they 
are  located  in  a  narrow  range  of  coplanar  distances  ( R ),  as  one 
would  expect  for  debris  from  a  dismpting  galaxy  in  a  nearly 
circular  orbit;  (2)  they  present  a  large  vertical  dispersion  of  ap¬ 
proximately  20  kpc  which  cannot  be  reconciled  with  a  disklike 
structure;  and  (3)  they  have  been  detected  in  positions  where  the 
fraction  of  thick-disk  stars  is  <  1 0%  of  the  total  number  of  stars  in 
that  direction.  The  only  doubtful  data  might  be  found  in  the  color- 
magnitude  detections  at  R  ~  13  kpc,  which  he  within  90%  frac¬ 
tion  of  thick-disk  stars.  These  points  correspond  to  what  our 
model  identifies  as  the  “close  ring,”  which  presents  a  higher  met- 
allicity  than  more  distant  detections. 

It  is  interesting  to  note  that  stars  stripped  from  the  Monoceros 
stream  progenitor  have  not  yet  contributed  to  the  thick-disk  pop¬ 
ulation.  The  decay  rate  of  the  dwarf  galaxy  appears  slower 
(around  3  kpc  in  3  Gyr)  than  the  mass-loss  rate  (50%  in  the  same 
time  period),  so  that  the  dwarf  galaxy  will  be  likely  destroyed 
before  reaching  the  inner  regions  of  the  thick  disc.  The  mass-loss 
process,  however,  depends  on  the  initial  structural  parameters  of 
the  dwarf,  which  were  fixed  ad  hoc  in  this  work.  More  accurate 
measurements  of  the  stream  surface  brightness  will  provide  better 
estimations  of  its  survival  time. 

6.  CONCFUSIONS 

In  this  paper  we  have  used  a  combined  semianalytic  TV-body 
technique  to  explore  the  nature  and  origin  of  the  Monoceros  tidal 
stream.  This  method  has  allowed  us  to  explore  a  large  parameter 
space  in  a  systematic  way,  with  the  goal  of  constraining  the  or¬ 
bital  properties  of  the  Monoceros  stream.  We  have  found  that  the 
available  observational  data  at  the  present  day  are  sufficient  to 
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robustly  constrain  some  aspects  of  the  Monoceros  stream  history 
and  progenitor.  In  particular,  (1)  the  heliocentric  distance  range 
at  which  the  stream  is  observed  plus  the  heliocentric  radial  ve¬ 
locity  curve  determines  an  orbital  eccentricity  ofe~0.10=L0.05, 
(2)  the  range  of  Galactic  latitudes  indicates  an  orbital  inclina¬ 
tion  of  i  ~  25°  =b  5°  from  the  disk  plane  and,  finally,  (3)  proper 
motions  of  stream  stars  are  only  compatible  with  a  prograde  orbit. 

We  found  that  it  is  considerably  less  straightforward  to  pre¬ 
dict  the  main-body  position  through  numerical  calculations.  The 
range  of  distances,  the  radial  velocity  curve,  and  the  vertical  ex¬ 
tension  of  debris  are  relatively  insensitive  to  the  final  location  of 
the  satellite  galaxy,  since  they  define  a  given  volume  in  phase- 
space  that  the  progenitor  fills  up  after  several  wraps.  Yet,  the  pro¬ 
jected  positions  of  stream  tails  are  time-dependent,  reflecting 
therefore  the  position  of  the  main  system,  or  equivalently,  the  final 
integration  time  of  our  TV-body  model. 

Besides  the  uncertainty  in  the  satellite  remnant  position,  the 
lack  of  detections  at  low  Galactic  latitudes  (owing  to  disk  absorp¬ 
tion  within  be  [—20°,  20°]),  as  well  as  for  a  complete  range  of 
Galactic  longitudes,  makes  it  difficult  to  constrain  the  location  of 
the  main  system.  Owing  to  these  limitations,  our  model  matches 
the  geometrical  and  kinematical  distribution  of  debris  if  the  main 
system  remnants  are  located  within  the  range  100°  >  l  >  200°, 
corresponding  to  integration  times  from  2.66  to  3.04  Gyr. 

On  the  other  hand,  observations  show  that  close  detections  of 
the  stream  have  higher  metallicities  than  the  distant  ones,  which 
provides  an  additional  constraint  to  our  model.  We  have  used  this 
constraint  by  assuming  that  the  progenitor’s  stellar  metallicities 
are  related  to  their  initial  binding  energy.  In  particular,  we  as¬ 
sume  that  stars  with  low  binding  energies  move,  on  average,  in 
the  outer  satellite  regions,  and  therefore  should  present  lower 
metallicities  than  those  with  high  binding  energies.  By  com¬ 
paring  the  initial  binding  energy  (metallicity)  distribution  with 
observations  for  the  range  of  locations  commented  above,  we 
find  that  the  metallicity  gradient  between  the  distant  and  the  close 
detections  can  only  be  reproduced  for  integration  times  longer 
than  2.9  Gyr.  In  that  case,  the  main  system  location  (/  ~  245°, 
b  ~  — 18°)  is  similar  to  that  of  the  CMa  dwarf,  with  a  helio¬ 
centric  distance  of  15.2  =L  0.8  kpc.  The  treatment  employed  to 
determine  the  distribution  of  metallicities  in  the  stream  is,  how¬ 
ever,  approximate,  because  we  have  no  way  of  knowing  the 
initial  metallicity  gradient  in  the  stream  progenitor.  Moreover, 
we  implicitly  assume  that  detections  with  different  metallicities 
reflect  an  initial  property  of  the  progenitor.  Other  scenarios  are 
possible,  and,  for  example,  we  cannot  rule  out  that  streams  with 
different  metallicities  at  different  distances  belong  to  different 
progenitors. 

The  best-fitting  model  predicts  a  halo  axis  ratio  of  qh  =0.6 
(oblate).  However,  the  poorly  constrained  selection  function  of 
observations  and  the  incomplete  area  coverage  leads  to  degen¬ 
erated  solutions  for  different  values  of  In  particular,  we  find 
that  similar  fits  can  be  obtained  with  halo  axis  ratios  0.6  <  qh  < 
0.8.  These  values  appear  in  agreement  with  the  result  of  Oiling  & 
Merrifield  (2000),  who  estimate  the  density  axis  ratio  of  the 
Milky  Way  to  be  around  0.8  from  kinematics  of  nearby  stars.  The 
Sagittarius  tidal  stream  has  also  been  used  to  constrain  the  shape 
of  the  Galaxy  halo.  In  that  case,  however,  the  results  are  con¬ 
tradictory,  and  so  whereas  Ibata  et  al.  (2001)  estimate  qh  ~  0.8  to 
—0.9  from  the  stream  width,  Helmi  (2004)  shows  that  the  recently 
obtained  radial  velocities  along  the  stream  can  only  be  explained  if 
the  Milky  Way  halo  is  fairly  prolate  (qh  —  1 .66).  Finally,  Johnston 
et  al.  (2005)  find  that  the  actual  Sagittarius  stream  precession  is 
only  consistent  with  an  oblate  Milky  Way  halo,  0.8  <  qh  <  0.9, 
ruling  out  the  prolate  case. 


We  must  emphasize  that  the  model  presented  here  is  far  from 
being  definitive.  Future  detections,  mapping  larger  areas  of  the 
sky,  will  provide  more  constraints  on  the  progenitor’s  orbit  and  a 
better  determination  of  parameters  listed  in  Table  3.  New  detec¬ 
tions  are  also  necessary  to  constrain  the  progenitor  location  by 
means  of  numerical  models. 

With  the  recently  available  kinematical  data  of  the  CMa 
dwarf,  which  Martin  et  al.  (2004a)  claim  to  be  the  progenitor  of 
the  Monoceros  stream,  we  have  integrated  the  orbit  backward 
in  time,  comparing  its  orbital  properties  against  those  of  our 
model.  We  find  contradictory  points  in  favor  of  and  against  that 
suggestion. 

In  favor: 

1.  The  orbital  eccentricity  of  CMa  (e  ~  0.16)  and  orbital 
sense  of  motion  (prograde)  are  consistent  with  those  of  our  model. 

2.  The  projected  location  of  CMa  (/  =  240°,  b  =  —8°)  is 
consistent  with  the  geometrical  and  metallicity  constraints  im¬ 
posed  by  the  Monoceros  stream  stars. 

3.  The  orbital  inclination  (/  ~  4^4  deg)  is  consistent  with 
that  of  our  model  within  1.5  cr  errors  from  proper  motions. 

4.  The  radial  velocity  of  CMa  (vr  =109  km  s-1)  is  also 
similar  to  that  of  our  model  (vr  =  135  km  s-1)  and  presents  a  low 
dispersion  incompatible  with  thick-disk  features. 

Against: 

1 .  The  orbits  of  globular  clusters  NGC  1851,1 904,  and  2298, 
which  Bellazzini  et  al.  (2004)  and  Martin  et  al.  (2004a)  claim  to 
be  an  evidence  that  the  CMa  overdensity  region  is  the  remnant  of 
the  Monoceros  stream  progenitor,  are  inconsistent  with  the  or¬ 
bit  of  CMa,  since  NGC  1851  and  NGC  1904  move  on  highly  ec¬ 
centric  orbits  and  NGC  2298  on  a  retrograde,  eccentric  one. 

2.  The  bimodality  in  the  velocity  distribution  of  CMa’s 
M-giant  stars  (Martin  et  al.  2004b)  can  be  reproduced  by  our  model 
in  regions  where  the  projected  positions  of  two  stream  tails  over¬ 
lap.  That  bimodality  in  vr  cannot  be  observed  in  the  central  part 
of  our  satellite  model  since  the  surface  density  is  several  orders 
of  magnitude  larger  than  that  of  overlapping  stream  tails.  This 
fact  seems  to  point  to  a  possible  contamination  of  background 
giant  stars  in  the  sample  used  to  measure  proper  motions  and  the 
radial  velocity  of  the  main  system.  As  a  result,  the  orbital  in¬ 
clination  might  have  been  underestimated. 

3.  The  main  body  of  our  best-fit  model’s  progenitor  is 
~  1 5  kpc  from  the  Sun — twice  as  distant  as  the  observed  CMa 
stellar  structure.  Since  the  mass  of  our  model  is  poorly  constrained, 
larger  mass  values  cannot  be  rejected.  In  particular,  dynamical 
friction  would  drive  our  satellite  model  down  to  7.2  kpc  in  3  Gyr 
for  an  initial  mass  of  approximately  1.8  x  109  M0,  whereas  our 
best-fitting  model  has  an  approximate  initial  mass  of  (6  =t  3)  x 
108  M0.  It  remains  unclear,  however,  whether  the  resulting 
distribution  of  debris  from  such  a  massive  satellite  undergoing 
tidal  disruption  would  also  reproduce  the  observations. 

Unfortunately,  with  the  results  obtained  in  this  work  we 
cannot  unambiguously  determine  whether  the  orbit  of  the  CMa 
dwarf  is  consistent  with  that  of  the  Monoceros  stream  progenitor, 
partly  because  of  the  uncertainty  in  the  available  observational 
data.  Stricter  criteria  to  select  stars  that  belong  to  the  CMa  dwarf 
appear  necessary  to  obtain  reliable  measurements  of  the  kine¬ 
matical  properties  of  this  system. 
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